************************************************************************
***************** Sanctuaries or Battlegrounds? ************************
************************************************************************

* Descriptive Statistics:
** Yearly
use "~\Sanctuaries or Battlegrounds Data - Yearly.dta"
estpost sum pro_gov_total rally mosque_center uni_students_pct public_jobs_pct NGO_total hospital_beds_per_1000 district_pop_total urban_pop_pct province_center_distance capital_distance
esttab using "~/Descriptive Stats.tex", label booktabs title(Summary Statistics\label {tab1})  cells("mean(fmt(2)) sd(fmt(2)) min max") replace noobs nomtitles nonum

** Monthly
use "~\Sanctuaries or Battlegrounds Data - Monthly.dta", replace
estpost sum  pro_gov_total rally mosque_center uni_students_pct public_jobs_pct social_protest anti_gov_protest sanction NGO_total hospital_beds_per_1000 district_pop_total urban_pop_pct province_center_distance capital_distance
esttab using "~/Descriptive Monthly Included.tex", label booktabs title(Summary Statistics\label {tab1})  cells(" mean(fmt(2)) sd(fmt(2)) min max") replace noobs nomtitles nonum


* Yearly models 
use "~\Sanctuaries or Battlegrounds Data - Yearly.dta", replace

* Create standardized version of variable of interests
foreach var in mosque_center uni_students_pct public_jobs_pct NGO_total hospital_beds_per_1000 district_pop_total urban_pop_pct province_center_distance capital_distance {
egen std_`var'= std(`var')
}
        

*Label the standardized variables
label variable std_mosque_center "Mosque centers (standardized)"
label variable std_uni_students_pct "University students (percent, standardized)"
label variable std_public_jobs_pct "State employees (percent, standardized)"
label variable std_NGO_total "NGOs (standardized)"
label variable std_hospital_beds_per_1000 "Hospital beds per 1000 (standardized)"
label variable std_district_pop_total "District total population (standardized)"
label variable std_urban_pop_pct "Urban population (standardized)"
label variable std_province_center_distance "Distance to province center (standardized)"
label variable std_capital_distance "Distance to capital (standardized)"


* Format the data for time-series analysis
xtset district_id_merge year

** First Model with Total number of pro-government mobilization events as the dependent variable
xtnbreg pro_gov_total std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance i.year,irr
estadd local year "Yes"
estadd local month "No"
estimates store model_1

** Second Model with number of pro-government rallies as the the dependent variable 
xtnbreg rally std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance i.year,irr
estadd local year "Yes"
estadd local month "No"
estimates store model_2


* Monthly model
use "~\Sanctuaries or Battlegrounds Data - Monthly.dta", replace


* Create standardized version of variable of interests
foreach var in mosque_center uni_students_pct public_jobs_pct NGO_total hospital_beds_per_1000 district_pop_total urban_pop_pct province_center_distance capital_distance social_protest anti_gov_protest {
egen std_`var'= std(`var')
}


* Label the standardized variables
label variable std_mosque_center "Mosque centers (standardized)"
label variable std_uni_students_pct "University students (percent, standardized)"
label variable std_public_jobs_pct "State employees (percent, standardized)"
label variable std_NGO_total "NGOs (standardized)"
label variable std_hospital_beds_per_1000 "Hospital beds per 1000 (standardized)"
label variable std_district_pop_total "District total population (standardized)"
label variable std_urban_pop_pct "Urban population (standardized)"
label variable std_province_center_distance "Distance to province center (standardized)"
label variable std_capital_distance "Distance to capital (standardized)"
label variable std_social_protest "Social protests (standardized)"
label variable std_anti_gov_protest "Anti-government protests (standardized)"


* Format the data for time-series analysis
xtset district_id_merge mdate

** Third Model with  monthly number of pro-government rallies as the the dependent variable 
xtnbreg rally std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance sanction std_social_protest std_anti_gov_protest i.month,irr
estadd local year "No"
estadd local month "Yes"
estimates store model_3

     

* Saving the models
esttab model_1 model_2 model_3 using "pro_gov_main_models_table.tex",  s(N year month aic bic , label("Num. Obs." "Year fixed effect" "Month fixed effect" "AIC" "BIC")) se cells(b(star fmt(2))  t(par fmt(2))) legend label replace eform noconstant star(* 0.05 ** 0.01 *** 0.001) drop(*year *month)

***-----------------------------------------------------------------------------------------------------------------------------------------***
* Models based on Actor Involvement
use "~\Sanctuaries or Battlegrounds Data - Yearly.dta", replace

* Create standardized version of variable of interests
foreach var in mosque_center uni_students_pct public_jobs_pct NGO_total hospital_beds_per_1000 district_pop_total urban_pop_pct province_center_distance capital_distance {
egen std_`var'= std(`var')
}
        

*Label the standardized variables
label variable std_mosque_center "Mosque centers (standardized)"
label variable std_uni_students_pct "University students (percent, standardized)"
label variable std_public_jobs_pct "State employees (percent, standardized)"
label variable std_NGO_total "NGOs (standardized)"
label variable std_hospital_beds_per_1000 "Hospital beds per 1000 (standardized)"
label variable std_district_pop_total "District total population (standardized)"
label variable std_urban_pop_pct "Urban population (standardized)"
label variable std_province_center_distance "Distance to province center (standardized)"
label variable std_capital_distance "Distance to capital (standardized)"

* Format the data for time-series analysis
xtset district_id_merge year

** Mosque Involvement Model
xtnbreg mosque_involved std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance i.year, irr
estadd local year "Yes"
estimates store mosque
** University Students Involvement Model
xtnbreg students_involved std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance i.year, irr
estadd local year "Yes"
estimates store students
** Public Employees Involvement Model
xtnbreg employee_involved std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance i.year, irr
estadd local year "Yes"
estimates store employees

* Actor Models Coefplots
grstyle init
grstyle color background white
grstyle gsize axis_title_gap tiny
coefplot mosque, keep (std_mosque_center) bylabel(Mosques)  ///
       || students, keep (std_uni_students_pct) bylabel(University Students) ///
	   || employees, keep (std_public_jobs_pct) bylabel(Public Employees)  ///
	   ||, xline(0) drop(_cons) byopts(row(1)) levels(95 50) ciopts(recast(. rcap)) ///
		   subtitle(, size(small) justification(center) ///
            bcolor(white) bmargin(small))
			
* Save the coefplot
graph export "~\Coefplot.png", replace


***-----------------------------------------------------------------------------------------------------------------------------------------***

*** Robustness checks:
		
* Per capita Transformation Models
* Yearly models
use "~\Sanctuaries or Battlegrounds Data - Yearly.dta", replace

** Create per capita variables and label them
gen mosque_center_per_capita=100000*mosque_center/district_pop_total
label variable mosque_center_per_capita "Mosque centers (per 100000)"
gen NGO_total_per_capita=100000* NGO_total/district_pop_total
label variable NGO_total_per_capita "NGOs (per 100000)"
gen district_pop_total_ln=ln(district_pop_total)
label variable district_pop_total_ln "District total population (log)"

* Format the data for time-series analysis
xtset district_id_merge year

** Yearly Models
xtnbreg pro_gov_total mosque_center_per_capita uni_students_pct public_jobs_pct NGO_total_per_capita hospital_beds_per_1000 district_pop_total_ln urban_pop_pct province_center_distance capital_distance i.year,irr
estadd local year "Yes"
estadd local month "No"
estimates store model_4

xtnbreg rally mosque_center_per_capita uni_students_pct public_jobs_pct NGO_total_per_capita hospital_beds_per_1000 district_pop_total_ln urban_pop_pct province_center_distance capital_distance i.year,irr
estadd local year "Yes"
estadd local month "No"
estimates store model_5

* Monthly Models
use "~\Sanctuaries or Battlegrounds Data - Monthly.dta", replace

** Create per capita variables and label them
gen mosque_center_per_capita=100000*mosque_center/district_pop_total
label variable mosque_center_per_capita "Mosque centers (per 100000)"
gen NGO_total_per_capita=100000* NGO_total/district_pop_total
label variable NGO_total_per_capita "NGOs (per 100000)"
gen district_pop_total_ln=ln(district_pop_total)
label variable district_pop_total_ln "District total population (log)"

* Format the data for time-series analysis
xtset district_id_merge mdate

xtnbreg rally mosque_center_per_capita uni_students_pct public_jobs_pct NGO_total_per_capita hospital_beds_per_1000 district_pop_total_ln urban_pop_pct province_center_distance capital_distance sanction social_protest anti_gov_protest i.month,irr iterate(50)
estadd local year "No"
estadd local month "Yes"
estimates store model_6

* Saving Per capita models
esttab model_4 model_5 model_6 using "pro_gov_appendix_pc_models.tex",  s(N year month aic bic , label("Num. Obs." "Year fixed effect" "Month fixed effect" "AIC" "BIC")) se cells(b(star fmt(2))  t(par fmt(2))) legend label replace eform noconstant star(* 0.05 ** 0.01 *** 0.001) drop(*year *month)

***-----------------------------------------------------------------------------------------------------------------------------------------***			
*Random-effects negative binomial models with robust standard errors 
* Yearly models 
use "~\Sanctuaries or Battlegrounds Data - Yearly.dta", replace

* Create standardized version of variable of interests
foreach var in mosque_center uni_students_pct public_jobs_pct NGO_total hospital_beds_per_1000 district_pop_total urban_pop_pct province_center_distance capital_distance {
egen std_`var'= std(`var')
}
        

*Label the standardized variables
label variable std_mosque_center "Mosque centers (standardized)"
label variable std_uni_students_pct "University students (percent, standardized)"
label variable std_public_jobs_pct "State employees (percent, standardized)"
label variable std_NGO_total "NGOs (standardized)"
label variable std_hospital_beds_per_1000 "Hospital beds per 1000 (standardized)"
label variable std_district_pop_total "District total population (standardized)"
label variable std_urban_pop_pct "Urban population (standardized)"
label variable std_province_center_distance "Distance to province center (standardized)"
label variable std_capital_distance "Distance to capital (standardized)"


* Format the data for time-series analysis
xtset district_id_merge year

** First Model with Total number of pro-government mobilization events as the dependent variable
xtnbreg pro_gov_total std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance i.year,vce(robust) pa irr
estadd local year "Yes"
estadd local month "No"
estimates store model_7

** Second Model with number of pro-government rallies as the the dependent variable 
xtnbreg rally std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance i.year,vce(robust) pa irr
estadd local year "Yes"
estadd local month "No"
estimates store model_8

* Monthly model
use "~\Sanctuaries or Battlegrounds Data - Monthly.dta", replace


* Create standardized version of variable of interests
foreach var in mosque_center uni_students_pct public_jobs_pct NGO_total hospital_beds_per_1000 district_pop_total urban_pop_pct province_center_distance capital_distance social_protest anti_gov_protest {
egen std_`var'= std(`var')
}


* Label the standardized variables
label variable std_mosque_center "Mosque centers (standardized)"
label variable std_uni_students_pct "University students (percent, standardized)"
label variable std_public_jobs_pct "State employees (percent, standardized)"
label variable std_NGO_total "NGOs (standardized)"
label variable std_hospital_beds_per_1000 "Hospital beds per 1000 (standardized)"
label variable std_district_pop_total "District total population (standardized)"
label variable std_urban_pop_pct "Urban population (standardized)"
label variable std_province_center_distance "Distance to province center (standardized)"
label variable std_capital_distance "Distance to capital (standardized)"
label variable std_social_protest "Social protests (standardized)"
label variable std_anti_gov_protest "Anti-government protests (standardized)"


* Format the data for time-series analysis
xtset district_id_merge mdate

** Third Model with  monthly number of pro-government rallies as the the dependent variable 
xtnbreg rally std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance sanction std_social_protest std_anti_gov_protest i.month,vce(robust) pa irr
estadd local year "No"
estadd local month "Yes"
estimates store model_9

* Saving Models
esttab model_7 model_8 model_9 using "pro_gov_with_robust_standard_errors.tex",  s(N year month aic bic , label("Num. Obs." "Year fixed effect" "Month fixed effect" "AIC" "BIC")) se cells(b(star fmt(2))  t(par fmt(2))) legend label replace eform noconstant star(* 0.05 ** 0.01 *** 0.001) drop(*year *month)

***-----------------------------------------------------------------------------------------------------------------------------------------***
* Model without Top Ten Populated Districts
* Yearly models 
use "~\Sanctuaries or Battlegrounds Data - Yearly.dta", replace

* Drop top 10 populated districts 
drop if district_en == "Tehran" | district_en =="Mashhad" | district_en =="Esfehan" | district_en =="Karaj" | district_en =="Shiraz"| district_en =="Tabriz"| district_en =="Qom"| district_en =="Ahvaz"| district_en =="Kermanshah"| district_en =="Urmia"

* Create standardized version of variable of interests
foreach var in mosque_center uni_students_pct public_jobs_pct NGO_total hospital_beds_per_1000 district_pop_total urban_pop_pct province_center_distance capital_distance {
egen std_`var'= std(`var')
}
        

*Label the standardized variables
label variable std_mosque_center "Mosque centers (standardized)"
label variable std_uni_students_pct "University students (percent, standardized)"
label variable std_public_jobs_pct "State employees (percent, standardized)"
label variable std_NGO_total "NGOs (standardized)"
label variable std_hospital_beds_per_1000 "Hospital beds per 1000 (standardized)"
label variable std_district_pop_total "District total population (standardized)"
label variable std_urban_pop_pct "Urban population (standardized)"
label variable std_province_center_distance "Distance to province center (standardized)"
label variable std_capital_distance "Distance to capital (standardized)"


* Format the data for time-series analysis
xtset district_id_merge year


* Yearly Model with Total Pro-Gov Mobilization as the Dependent Variable
xtnbreg pro_gov_total std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance i.year,irr
estadd local year "Yes"
estadd local month "No"
estimates store model_10

** Yearly Model with number of pro-government rallies as the the dependent variable 
xtnbreg rally std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance i.year,irr
estadd local year "Yes"
estadd local month "No"
estimates store model_11

* Monthly model
use "~\Sanctuaries or Battlegrounds Data - Monthly.dta", replace

* Drop top 10 populated districts 
drop if district_en == "Tehran" | district_en =="Mashhad" | district_en =="Esfehan" | district_en =="Karaj" | district_en =="Shiraz"| district_en =="Tabriz"| district_en =="Qom"| district_en =="Ahvaz"| district_en =="Kermanshah"| district_en =="Urmia"

* Create standardized version of variable of interests
foreach var in mosque_center uni_students_pct public_jobs_pct NGO_total hospital_beds_per_1000 district_pop_total urban_pop_pct province_center_distance capital_distance social_protest anti_gov_protest {
egen std_`var'= std(`var')
}


* Label the standardized variables
label variable std_mosque_center "Mosque centers (standardized)"
label variable std_uni_students_pct "University students (percent, standardized)"
label variable std_public_jobs_pct "State employees (percent, standardized)"
label variable std_NGO_total "NGOs (standardized)"
label variable std_hospital_beds_per_1000 "Hospital beds per 1000 (standardized)"
label variable std_district_pop_total "District total population (standardized)"
label variable std_urban_pop_pct "Urban population (standardized)"
label variable std_province_center_distance "Distance to province center (standardized)"
label variable std_capital_distance "Distance to capital (standardized)"
label variable std_social_protest "Social protests (standardized)"
label variable std_anti_gov_protest "Anti-government protests (standardized)"


* Format the data for time-series analysis
xtset district_id_merge mdate

** Monthly Model with number of pro-government rallies as the the dependent variable 
xtnbreg rally std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance sanction std_social_protest std_anti_gov_protest i.month,irr
estadd local year "No"
estadd local month "Yes"
estimates store model_12

* Saving Models
esttab model_10 model_11 model_12 using "pro_gov_without_populated_districts.tex",   s(N year month aic bic , label("Num. Obs." "Year fixed effect" "Month fixed effect" "AIC" "BIC")) se cells(b(star fmt(2))  t(par fmt(2))) legend label replace eform noconstant star(* 0.05 ** 0.01 *** 0.001) drop(*year *month)

***-----------------------------------------------------------------------------------------------------------------------------------------***

* Linear Regression Model for Raisi's Vote Prediction
use "~\Sanctuaries or Battlegrounds Data - Yearly.dta", replace

* Create standardized version of variable of interests
foreach var in mosque_center uni_students_pct public_jobs_pct NGO_total hospital_beds_per_1000 district_pop_total urban_pop_pct province_center_distance capital_distance {
egen std_`var'= std(`var')
}
        

*Label the standardized variables
label variable std_mosque_center "Mosque centers (standardized)"
label variable std_uni_students_pct "University students (percent, standardized)"
label variable std_public_jobs_pct "State employees (percent, standardized)"
label variable std_NGO_total "NGOs (standardized)"
label variable std_hospital_beds_per_1000 "Hospital beds per 1000 (standardized)"
label variable std_district_pop_total "District total population (standardized)"
label variable std_urban_pop_pct "Urban population (standardized)"
label variable std_province_center_distance "Distance to province center (standardized)"
label variable std_capital_distance "Distance to capital (standardized)"


* Format the data for time-series analysis
xtset district_id_merge year

* Linear regression model
reg Raisi_votes_2017_pct std_mosque_center std_uni_students_pct std_public_jobs_pct std_NGO_total std_hospital_beds_per_1000 std_district_pop_total std_urban_pop_pct std_province_center_distance std_capital_distance if year==2016
estimates store reg_model

* Saving Linear Regression Model
esttab reg_model using "pro_gov_reg_model.tex", aic bic se cells(b(star fmt(2))  t(par fmt(2))) legend star(* 0.05 ** 0.01 *** 0.001) label replace

